STABILITY ANALYSIS OF FLOCK AND MILL RINGS FOR 2ND 
ORDER MODELS IN SWARMING 
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^^ Abstract. We study the linear stability of flock and mill ring solutions of two individual 

Cn based models for biological swarming. The individuals interact via a nonlocal interaction 

^ potential that is repulsive in the short range and attractive in the long range. We relate 

Mh the instability of the flock rings with the instability of the ring solution of the first order 

-^T model. We observe that repulsive-attractive interactions lead to new configurations for the 

flock rings such as clustering and fattening formation. Finally, we numerically explore mill 

patterns arising from this kind of interactions together with the asymptotic speed of the 

system. 



OS 



cd 



1. Introduction 



CM 
< 

4^ Individual-based models (IBMs) appear in biology, mathematics, physics, and engineering. 

^ They describe the motion of a collection of N individual entities, so the system is defined 

^ on a microscopic scale. IBMs are good models for some applications when the number of 

I— ' particles is reasonable. Nonetheless, when the number is large, it is more reasonable to use a 

continuum model. Some continuum models, like the one in [131 US], are derived as a mean- 

^ field particle limit leading to a mesoscopic kinetic description of the problem. At this level, 

Q\ one looks at the probability density of finding particles at a certain position and velocity at a 

^ given time. Several models have been proposed to describe the flocking of birds [HI EHl (H [27], 

ly^ the formation of ant trails [18], the schooling of fish [2Tl[7l|5], swarms of bacteria [23], etc. 

I* These models can include some rules that describe the behavior of each individual of the 

^^ system. Such mechanisms can help to describe the influence of each individual on the others, 

CO depending on their relative position and velocity. For instance, one example is the classical 

^^ three zone model [H [22] . A three zone model describes how social the individual is in the 

^ following sense. If two individuals are too close, they want to have their own space (repulsion). 

S>( When one individual is far from the group, it wants to go back and socialize (attraction). And 

%^ finally, in the group, each individual tries to mimic the behavior of the others (orientation). 

Other models just consider rules for orientation, like the Vicsek model [291I16]. In this case, 
there is a mechanism of self-propulsion in which each individual moves with constant speed 
and adopts the average direction among their local neighbors. 

We focus our study in the analysis of two particular examples of IBMs. The first one is a 
self-propelled interacting particle model that was introduced in [26] and extensively studied 
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We are going to consider the same self-propulsion/ friction term used in [17^ .13]. 

S{\vj\) = a- P\vj\^, a,/3>0. 

Note that such a term gives us an asymptotic speed for the particles, equal to ^ aj ^. In 
these references, the authors study ([l]) with pairwise interaction given by the so-called Morse 
potential 

with Cai Cr denoting the attractive and repulsive strengths and Ia^ Ir their respective 
length scales. They find and describe several patterns for the asymptotic behavior in 2D. 
They observed flocking behavior, mill on a ring, and clustering when particles are milling. In 
[9], a well-posedness theory is developed for ([l]) proving the mean- field limit under smoothness 
assumptions on the potential. The authors show convergence of the particle model toward a 
measure solution of the kinetic equation. 

We perform an analysis on the stability of flock rings and the mill rings as asymptotic 
solutions for Q. The ring solution was recently studied in [25 ^ 6J where the authors in this 
work do a careful general linear analysis of the rings for the first order model 
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Xj = Y,^W{X,-Xi), j = l. 



.N. 
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Part of the analysis of ^ is used to study the the stability of mill rings in ([T]), described 
in [T7||. Related pattern formation in the associated first order model has been studied in 

mm- 

Another second order model that we are going to study is 



N 



. N , ., 

- 2 H{x, - xi){vi - V,) + - 2 ^W{xi 
1 



N 






j = h...,N 



(3) 
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with Xj,Vj G R^ where the velocity Vj is described by the Cucker-Smale alignment term H 
and the pairwise interaction by a repulsive- attractive radial potential W{x) = k{\x\). 

Even if the analysis has been done in full generality for the parameter functions of the 
model H and VF, we will emphasize the results in some relevant cases. For instance, we 
consider the case of power law repulsive-attractive potentials [251 El 



k{r) = 



a > b > 0. 



(4) 



For the Cucker-Smale alignment [111 HSl [201 [13 HO] , ^ relevant case is H{x) = fifda^l) with 

1 



gir) 



(l+r2)7' 



7>0. 
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The main results of this work show that the flock ring is unstable in the second order 
models ^ and ([s]) if and only if its spatial shape is unstable in the first order model Q. We 
use the same kind of strategy as in [6]. 

The main idea is to study the stability of the system of ODEs ^ by analyzing the eigenval- 
ues of a suitable linearization with restricted perturbations. We consider particular perturba- 
tions of the flock rings in such a way that translational invariance is avoided while preserving 
the mean velocity. This is really needed since the linearized system associated to ([l]) is al- 
ways linearly unstable due to translations. Translational invariance implies the existence of a 
generalized eigenvector associated to the zero eigenvalue of the matrix defining the linearized 
system. We characterize all cases in which the linearized system has eigenvalues with zero 
real part and their consequences in the instability condition. An analysis of the stability of 
the family of flock solutions is under way in [llj. 

In addition to flock rings, other spatial shapes are possible as asymptotic solutions. One 
can also observe flocks on annuli (fattening), lines or points (clustering). These patterns can 
be explained due to the results in pj. 

For the mill ring analysis, we start from the results of [6J to explore other mill configurations 
that appear with repulsive-attractive potentials. We numerically investigate the formation of 
fat mills due to the repulsive force and the formation of clusters when varying the asymptotic 
speed. In addition, we show some switching behaviors between flock and mill rings. 

The structure of the paper is as follows. In Section [2] we study the microscopic and the 
mesoscopic models and we give the definitions of the main objects in our studies, the flock 
and mill rings. In Section |3] we do a linear stability analysis on the flock rings for models 
^ and ([3]). We also explore the fattening and cluster formation. Finally, in Section [4] we 
do a similar study on the stability for mill rings. In all sections we have performed several 
numerical tests supporting our theoretical results. 

2. MICROSCOPIC & Mesoscopic models 

Let us introduce some particular solutions of the particle model ([l]) and its continuum 
counterpart. 

2.1. Flock and mill solutions: microscopic model. 

Definition 1. We call a flock ring^ the solution of ^ such that {xj}^^-^ are equally distributed 
on a circle with a certain radius, R and {vjW^-^ = uq, with \uo\ = ^ aj ^. 



Definition 2. We call a mill ring^ the solution of ([!]) such that {xj}^^^ are equally distributed 

on a circle with a certain radius, R and {vj}^^^ = u^ = -yJa//3 x^/\xj\ with x^ the orthogonal 
vector. 



By abuse of notation, we will write \uo\ for u^ since u^ = \J olJ ^ for all j = 1, . . . , A^. 

Moreover, we will make use of notation \u^\ for both flock and mill rings indistinctly. 

All over the paper, we will identify e^^ = (cos 0, sin 0) and use x to identify the vector and 
the complex numbers indistinctly 

x,{t) = R (cos (^i + ^i) ,sin {^^j^^i\\ = i?e^^e-*. (5) 
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Figure 1. Flock and mill ring solutions. 
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In the case of mill rings, we are looking for a solution of the form ([s]). The case of a flock 
ring in the comoving frame is equivalent to looking for a solution of the form ^ with cj = 0. 
Plugging ([5]) into Q, we obtain 



N 



Y, VH^(x, - xi) = 0; 



1^1 
1^3 



therefore, R is determined only by the repulsive-attractive potential. 
In the case of mill rings, we have 



27rj 
Vj\i) = Xj\i) = nujie ■" e , ^j = 

and thus, R'^oj'^ = ^. Moreover, by taking the derivative 



Vj{t)=xjit) = Rioie''^e'^\ 0j = ^, 



1 ^ 



-oj'^Xi = —up'—- y '(x/ — Xn) = 0. 






Plugging this into ([T]), we get 



N N 

J] [VW{xi - xj) - Lj\xi - xj)] = J] VW{xi - xj) = 0, 
1=1 1^1 

with W{x) = W{x) — cj^^. Thus in order to find the radius for flock and mill rings, we need 
to solve the same equation. This expression implies that the spatial shape has to balance 
attraction versus repulsion and centrifugal forces. Now, a direct computation yields 

\xj — xi\ = 2i?sin f — 

for all times. One can easily compute that 

sin(^) cos(^)\ (cos {9 



xi-Xj = 2Rsm(^ \ I I I , p = I - j, 



cos(f) sin(f)y \sm{9j 
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and 
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N 



Y.'^Wilxj-xiD^Y^ia 



Xl 



k'{\xj-xi\) 



1=1 

1^3 



I '^3 '^l I 

sin (2^) cos (Oj) + cos (f ) sin {Oj] 



1=1 

N 

1=1 V cos (f ) cos (dj) + sin (f ) sin (9,, 

By symmetry we can assume j = N, and thus 

AT N-1 /-sin [^y 



\k\\xj -xi\). 



J]VW{\x,-xi\)= J] 



1^1 

1^3 



COS 



\N ) 



\kWxj -xi\), 



changing j to N — j, we finally obtain 



N 



^VWi\xj-xi\) 



p=i 



1=1 

1^3 
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(6) 
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As a conclusion, the radius of a flock or mill ring solution is characterized by 



p=i 



\nJJ 



x^n^C^)-^^ (7) 



For general potentials there can be more than one flock or mill solution. In the case of the 
power law potentials Q, there is only one solution. Condition ^ reads 

p—0 p=0 p=0 

To prove uniqueness, we notice that the function /(r) = Cir^ — C2r^ — C3 with a > 6 > 
and Ci > C2 > 0, C3 > 0, has only one zero. Computing the first derivative and looking for 
critical points, we obtain ri = and r2~ = ^^. Taking the second derivative and evaluating 

at r2, one obtains f"{r2) = ^2" C2b{a — 6) > 0, so r2 is a local minimum. For < r < r2, one 
has /^(r) < whereas for all r e (r2, +oc), one has /'(r) > 0. Then we conclude that /(r) 
has a unique zero. Notice that the solution to ^ depends on the number of particles and we 
will use the notation R = R{N). 

2.2. Flock and mill solutions: mesoscopic model. In this subsection we characterize 
the radius of flock and mill rings. We introduce the function 



i^a{s) 



1 r 

^ Jo 



(1 - 5cos6>)(l + 5^ - 25cos6>)^ d0, 



already analyzed in [3] . A change of variables in the previous function shows that 



V'a(l) 



-)a—l 



-B 



TT 



a + 1 1 
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Lemma 1. If N ^^ co then R{N) -^ ^a^d'^ol) where i^^^d'^ol) ^^ the solution of the following 
equation: 

Proof We first take 2^"-*^-^ V Jq sin^(^). Multiplying and dividing by tt we obtain the 
following equality 

N-l 






v-1 /r)7r\\ 1 r^ 1 / r^/2 \ 

J] sin« [j-j = 2^-1- sin^(x) dx = 2^'^- 2 sin^(x) dx . (9) 



Now, we use the following expression for the Beta function 

r7T/2 



rTTfZ 

Jo 
with X = ^^ y = ^^, and using that B{x^y) = B{y^x) in ^ together with ([s]) to obtain 

lim 2-iif^ysin«f^)V2^-'-^f^^.-l =^.(1)- 
AT-oo ^l^4f, \N)) tt V 2 '2y ^""^ ^ 



p=0 

The same reasoning works by changing a for h in the second term in Q. For the third term 
we use the fact that we can compute the exact sum 

I S -n^ m = 1. 

n 

Remark 1. In the case of flock rings u = 0, their radius is determined by the radius of the 
aggregation ring found in ^ 



1 B{'-^,1) ^ "- 

2 1,5(^,1) 



R{N) ^ Rab = - \ -^-Ij^{ ] as N^^. 



Remark 2. Let W{x) = k{\x\) be a general interaction potential. Call f{r) = —k'{r)/r. 
Then the radius of the ring is determined by 

' /(2i?sin(0))sin2(0)d^ = O, 
Jo 

as shown in [6J. 

Remark 3. The corresponding mesoscopic model to the particle system ^, as proven in ^, 
is given by the kinetic equation 

j^ + v-V:cf + div4(a - I3\v\^)vf)] - div^[(Va.VF * p)f] = 0, (10) 

where 

P{t,x) = f{t,x,v)dv. 

Jr2 

It was shown in [12j that singular solutions of the type 

(id X \ 
\ P X J 
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with p{t^x) the uniform distribution on a ring, are weak solutions of the kinetic model (10)^ 
called the flock and mill ring continuous solutions respectively. 



3. Linear stability analysis for flock rings 

We will now focus on the stability analysis of flock rings for some particular perturbations 
in terms of the parameters of the model (a, 6, u^). We take advantage of the careful stability 
analysis of the ring solutions of the aggregation equation performed in [6J. 

3.1. Stability of flock solutions without the Cucker-Smale term. We consider the 
model ([l]) and we perform the change of variables to the comoving frame 



j = l,...,Ar, 



\yj{t) = Xj{t) -uot 

\zj{t) =Vj{t) -uo 

where uq is the asymptotic velocity of a fixed flock ring. Therefore the system ([l]) reads 

r d 



(11) 



dt 
{ d^ 

dt 



Vj = ^3 - «o 



N 



Zj = {a- I3\zj + uo\^){zj + ^^o) - ^ XI ^^(^i ~ Vk) 



l,...,N. 



SoQzjl) 



k=l 
kitj 



A flock ring can then be characterized as a stationary solution of the form iy^,z^ 
(i?e*^J, O) , where 6j = '^ ioi j = 1, . . . , N . This stationary solution satisfies 

5o(|4|) = 0, V5o(|4|) = -2/3nono*, = ^ ^'(bi " 2/°!)^^" ^'^ 



As in [6] , we restrict the set of possible perturbations of the flock solution to those of the form 

y,it) = Re'''il + h,it)), 



where hj e C, such that \hj\ « 1 and 



N N 

j=i i=i 



(12) 



The first restriction is to avoid the zero eigenvalue due to translations. The second one comes 
from the fact that the mean velocity of the perturbed system should be uq. More general 
perturbations will generically lead to other flock solutions with different asymptotic velocity 
Uq. Their orbital stability will be analyzed elsewhere [11]. Therefore, the perturbed system 
reads 

^j^y,{t) = Re^'^h'^ = ~z, 



1=1 
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The linearization of the system around the flock solution {y^, z^) reads as 



^^J,{t) = Re^'^h'^=zj 



dV 



1 



N 



Re^'^h] = -2/3non* - j^J] VW{yj - m) 



.N. 



From the previous equation, we can characterize h'^ as 

h'; = 2 [GiW2)ihj - e'%) + G2i<i>/2){hi - e'%)]\ - 2(3uoulh'^, 



where 



1^3 

^ j^ and 



Gi{4>) =^ [-a{2R\ sin0|)«-2 + h{2R\ sin^l)''-^] , 

G'2(</') =^ [-(a - 2)(2R\ sin^|)«-2 + (6 - 2){2R\ sin0|)''- 

for the power law potentials. The details of these previous computations for general potentials 
can be found in [6]. 

Let us consider the following ansatz for perturbations hj 



hj=^+it)e^^''+Ut)e 



-imOj 



m = 2,3, . . , 



(13) 



which satisfles conditions ( |12[ ). We need to exclude the case m = 1 since it leads to a zero 
eigenvalue due to the rotational invariance of the system. Following the same strategy as in 
|fij and some computations, we flnally deduce 



C ~ [him) hi-m) U- 



2/3uoUq 



T ii+ 



M 



with /i and I2 real functions given by 



NI2 



1^3 



N/2 



him) =J]G2 i<t>m ie'^^ - e^^) = 4 J] G2 (^ 






(14) 
(15) 



The previous system can be written also in the following form 



d 
dt 






Id \ 

M -2(3uou^) 



/CA 




(U\ 


k- 
v+ 


= L 


i- 
v+ 



(16) 



where (r?+,r?_) = i^'+,0. 
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If we do not assume (12) and (13), then we cannot reduce the analysis to a 4 x 4 system. 



An arbitrary perturbation for general flocks leads instead to a matrix of the form 




where the partition into 2N x 2N sub-blocks reflects the distinction between position and 
velocity contributions to the Jacobian: The symmetric matrix M is the 2N x 2N Hessian 
that results from linearizing the first order system ^ about a given flocking configuration, 
whereas U denotes a block-diagonal matrix with N blocks of the 2x2 matrix i^o'^o^ along the 
diagonal. By rotational invariance we can reduce to the case U{) = ei = (1,0), so that the 
block matrix U acts on x = (xi, . . . , xn)^ g R^^, xi e R^, according to the relation 

(t/x), = (^^''^Q^l' 

We now turn to the task of characterizing the eigenvalues of L in terms of the eigenvalues 
of M. In other words, we aim to characterize the stability of a flock in terms of the stability 
of its spatial shape as a solution to the first order model. To fix the notation, we write the 
eigenvalue problem for the flock as 

Af!i=(:. jiV!)=i^(j). (17) 



vy) \m -2pu) \W 



where the matrix M determines the stability of the flocking configuration as a solution of the 
first order model. For any given eigenvector (x, v) e C^^ x C^^ of the full system (17), we 



always assume the normalization x*x = 1. Substituting the first equation Ax = v into the 
second equation yields the equivalent statement 

A^x + 2/3At/x - Mx = 0. (18) 

Let |x|2 denote the semi- norm on C^^ defined according to 

N 



and let E^ ^ C^ denote the subspace 



i=l 



E^ := {x G C^^ : |x|2 = 0} = ker(t/). 

Premultiplying by x*, the fact that x*t/x = |x|2, the normalization on x and the quadratic 
formula combine to imply the key identity 



A = -/3|x|^ ± Y/32|x|| + x*Mx. (19) 

As M is symmetric, we may write its 2N real eigenvalues and corresponding normalized 
(x*x = 1) eigenvectors as 

l^2N ^ l^2N-l ^ • • • ^ M2 ^ Ml Mx^ = /X^Xi. 

The notation aL(A), aM(M) will denote the algebraic multiplicities of A, /x as eigenvalues of 
their respective matrices. The bulk of the analysis lies in characterizing the eigenvalues A of 
the full system (fl7]) that have 3ft(A) = 0. 
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Lemma 2. Let A denote an eigenvalue of (17). Then 3ft(A) = and Q'(A) / z/ and only 
if X = ±.i^J—|Jik for some k with j^k < and ^k ^ E^ . The eigenspace consists only of 
eigenvectors. 



Proof If :sik e E^^ then (18) reads A^x/^ = Mx/^, or equivalently A^ = /i^^. To have Q^(A) / 



then requires /x/e < 0. Conversely, if (x, Ax) denotes an eigenvector with $ft(A) = and 



^{X) / 0, the formula (19) implies that necessarily x e E^ ^ and therefore Mx = A^x. Thus 
A^ = fik for some /ik < 0. 

To show the last statement, suppose a generalized eigenvector existed that is not an eigen- 
vector. Then there exists an eigenvector (x. Ax) with x e E^ so that the system of equations 

-Aid Id A fu\ /x ^ . . 



M -2I3U - XldJ \wj \X^^ 

has a non-trivial solution. Substituting the first equation w = Au+x into the second equation, 
then pre-multiplying by x* demonstrates 

Mu - 2/3t/w = 2Ax + A^u 

x*Mu = 2A + A^x^u. 

The last line follows as x*x = 1 and x e E^ = ker(t/). The symmetry of M and the fact that 
Mx = A^x combine to show x*Mu = A^x*u. Thus A = 0, leading to a contradiction. D 



Lemma 3. Let /3 > 0. Then X = is an eigenvalue of ( |17| ) and (x, 0) is a corresponding 
eigenvector if and only z/ Mx = 0. If ^ ^ E^ then (x, 0) generates a single generalized 
eigenvector, whereas if^^ E^ then (x, 0) generates no generalized eigenvectors. 



Proof. The first statement follows trivially from (18). To see the second statement, consider 



the system of equations (20) with A = 0. This reduces to the equations w = x and 

Mu = 2/3t/x, 

which by premultiplying by x* as before and using the fact that Mx = necessitates x e E^ 
as /3 > 0. If indeed x e E^ then any u e ker(M) suffices. Without loss of generality, take 
u = X itself. If (x, 0) generates a second generalized eigenvector then the system of equations 

Id \ /u\ ^ /x' 

^M -2PU) \^) ~ Vx^ 

has a non-trivial solution. As then w = x and x e E^ this reads Mu = x. Premultiplying 
one last time by x, the facts that Mx = and x*x = 1 combine to produce the contradiction 
= 1. D 

This lemma yields, as a corollary, the algebraic multiplicity aL(0) of zero as an eigenvalue 
of the second order system. 

Corrolary 1. Let /3 > 0. Then 

aL(0) = dim(ker(M) n E^) + dim(ker(M)). 

Let aM,i(0) := dim(ker(M) n E^)^ so that aL(0) = aM,i(0) + aM(0). Note that neither 
quantity depends on /3, and the conclusion holds whenever j3 is positive. Thus, if /3 e (0, oo) 
it follows that aL(0) is constant. Moreover, Lemma [2] holds uniformly in j3 as well. Let 
ii < i2 < ' • ' < ik ^ 2A^ denote those (possibly non-existent) indices where fii^ < has an 
eigenvector x^. e E^ . The two lemmas then combine to show: 
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Corrolary 2. Let /3 > 0. Then 

det(L - Aid) = x<^^M0HaM0)iik^^^x^ _ ^^^)p^(X), 

The roots of the polynomial pi^{X) all have non-zero real part. 



This corollary, along with the formula ( |19[ ), suffice to estabhsh the desired result: 

Theorem 1. The linearized second order system around the flock ring solution ^ has an 
eigenvalue with positive real part if and only if the linearized first order system around the 
ring solution has a positive eigenvalue. As a consequence, the flock ring solution is unstable 
for m-mode perturbations for the second order model ^ if and only if the ring solution is 
unstable for m-mode perturbations for the first order model Q. 

Proof. Suppose first that /xi ^ 0. Then x*Mx ^ for any x, whence all eigenvalues A of L 
have non-positive real part due to (fT9|). Conversely suppose /xi > and let A denote the set 



A:= <I3e[0,oo): max 3ft(A) > 
Xe(j{L) 



Note that e ^ due to (19). Indeed, then (xi, y^^xi) defines an eigenvector with eigenvalue 
^ = v^ ^ ^- ^y continuous dependence of the eigenvalues of L on ^, it follows that A is 
relatively open. To show that it is also relatively closed, let I3i e A and /3/ ^- /3o e (0,oo). 
Up to extraction of subsequences, it follows that there exists a corresponding sequence A^ of 
eigenvalues with 3ft(A^) > converging to some Aq with 3ft(Ao) ^ 0. Moreover, by continuous 
dependence of the coefficients of ^/^(A) on /3, the roots of p/3^(A) converge to roots of P/3q(A). 
Thus P/3q(Ao) = 0. As no such root can have zero real part by corollary l2| 3ft(Ao) > and 
/3o e A. As .4 / it follows that A = [0, oo) as desired. The last part of the theorem is a 
direct application of the first part to the 4x4 Tn-mode perturbation matrix in (|16|) . D 



Remark 4. As an artifact of translation invariance in the first order model, the vector defined 
by e2 := (0, 1, ... , 0, 1)-^ g M?^ always defines an eigenvector o/M with eigenvalue zero. Due 
to the fact that e2 e £'^, Lemma [s] implies that (e2,e2) furnishes a generalized eigenvector 
with eigenvalue zero, so that the flock is always linearly unstable for the model ^ . 

3.2. Numerical validations. In this section, we perform some numerical computations to 
show stability regions for the flock ring. Moreover, we will show the formation of clusters and 
the fattening instability. Due to Theorem [T] and Corollary [4j applied to the 4x4 matrix in 



(16), we are reduced to study the determinant of the matrix M for the flock stability. Note 
that for fixed values of N and m the determinant of M is a function of the parameters a and 
b and such that 

D{a,b) :=det{M) = /i(m)/i(-m) - (hMf^ 
T{a,b) :=trace(M) = /i(m) + /i(-m). 

Remark 5. Using the results of [6, Theorem 3.1] one is able to estimate the asymptotic value 
of the determinant of M. In our case, using W{x) = ^ ^ one obtains that 

det(M) - Cm~^^^ as m^ oo, 

where C > and be (1, 2) u (4, 6) u (8, 10) u • • • . In these cases det(M) > and trace(M) < 0. 
Moreover, this result shows that there is no spectral gap for large modes m since det(M) -^ 

as 777. ^- 00. 
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In Figure [2] we recover some results on the stability already shown in [25j . Since Ii{m) 
and l2{Tn) depend on the powers a and b of the power law potential, we plot in the parameter 
region {(a, 6) : a > 6 > 0} the stability and instability regions depending on the determinant 
of M and its trace. We show the cases tti = 3, 4, 5 for a fixed N = 100000. In Figure |3] we 
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Figure 2. Stability regions for different perturbations modes m = 3,4,5. 
Stability regions correspond to a couple (a, 6) of parameters where D{a,b) > 
0, r(a,6)<0. 

compute the stability area as a function of a and b. To do so, we compute the intersection 
of all stability areas for m ^ 2. It can be observed from our tests that the stability area 
shrinks when the number of particles increases. Moreover, it is observed that in the limit 
when N ^ GO, the lower boundary of the stability region converges to the dashed line. The 
red dashed line is the curve b = ^^ that corresponds to the tti = +00 mode. This curve 
is the separatrix of the ins/stability regions for the continuous delta ring of the first order 
continuum model, studied in [25113]. 
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Figure 3. Stability areas for flock ring solutions for different values of N . 
From left to right: N = 100,1000,10000. Markers {;^) and (0) indicate the 
explored parameters respectively in Table [l] and Table [2j 



3.2.1. Cluster formation. The formation of clusters occurs when the repulsion strength is 
small. In other words, this phenomenon depends on how singular the potential is at the 
origin. We show the bifurcation diagram for the phase transition between equally distributed 
flock and flock with cluster formation. Figure |4] is performed using N = 1000 particles 
equally distributed on the stable circle with all the velocities aligned. We let them evolve 
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Figure 4. Bifurcation diagram for cluster formation at Tf = 500, with N 
1000 particles, a = 5, |i^o| = 2.5. 



until Tf = 500. We fix the parameters \uq\ = 2.5, a = 5 and vary 6 along the axis. The 
vertical axis represents the increment of the relative errors 

u Nil _ llM - Mo l|2 

Up llre^ — II ^11 

with increasing 6, where jj^q^ is the uniform distribution along the stable ring of N particles 
and /i^ the distribution at time Tf. Simulations are performed with MATLAB and the 
evolution of the system of odes is solved with the ode^B routine with adaptive time step. 
Table [l] illustrates different possible final states for different choices of the parameters a, h. 



a = 3, 6 = 2.5 



a = 5,6 = 4.1 



a = 7,6= 1.5 





^^^^^ 




Table 1. Long time simulations with A^ = 1000 particles. The location of 
parameter values are marked as {^) points in the central plot of Figure |3J 



3.2.2. Fattening formation. We show the transition diagram between a flock on a ring and a 
flock on an annulus. In this case, the fattening phenomenon occurs when the parameters of 
the potential cross the lower boundary of the stability region. We numerically characterize 
this behavior in a similar way as in the previous subsection. 
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Figure 5. Bifurcation diagram for fattening instability at Tf = 500 with 
N=1000 particles, a = 5, |i/o| = 2.5. 



Figure [5] is performed using N = 1000 particles equally distributed on the stable ring 
already in the steady state with all the velocities aligned. We then let them evolve until 
Tf = 500. We fix the parameters \uq\ = 2.5, a = 5 and vary b along the axis. The vertical 
axis represents the increment of the relative errors 



ll^^llre^ 



II A^ A^ll 



H 



N\\ 



for increasing 6, where 770^ represents the average distance from the center of mass for N 
particles in a flock ring formation, i.e., 77q^ = R and 77^ is the average distance from the 
center of mass at time Tf. Table [2] shows the final states for some particular choices of 
parameters after stabilization. 



5,6= 1.1 



5,6 = 0.5 



7,6 = 0.5 






Table 2. Long time simulations with N = 1000 particles. The location of 
parameter values are marked as (0) points in the central plot of Figure |3J 
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3.3. Stability of flock solutions with the Cucker-Smale alignment term. As in Sec- 
tion 3.1 we perform a study on the flock stability for model Q. If we use the same change of 



variables as in (11), then system ^ reads 



Vj = ^j - ^0 = Zj 

N 



N N 

^J = ]^J]H{\yi-y,\){zi-Zj)^-J]VW{yi 



Vj] 



j = l,...,7V (21) 



-1 

1^3 



We give a characterization of the flock solution in the complex plane for (21) with (^^, z^-) = 
(i?e^^-^, 0), where 9j = -jf. We consider then the perturbed solution 

with hj such that \hj\ « 1 and satisfying ( [T2| ). Considering the following relations 
yi-y^=Re'^^ (^e^^hi-hj^ 



\yi -yj\- 2i? 



sm 



R 



sm 



1 - e^*-) {hi - hj) + (l - e-^*-) {hi - hj) 



zi - Zj = Re'^^ (e^Ph'i - h'j) , 



where (pp = 2n{l — j)/N = 2'Kp/N. We linearize the Cucker-Smale alignment term around 
the solution up to first order, leading to 



H{\yi-yj\)^H{2R\sm{<Pp/2)\) 

+ H'{2R\sm{<l)p/2)\)- 



R 



sin (%) 



1 - e^"^^) {hi - hj) + (l - e-''t>^^ {hi - hj) 



Substituting the linearization in (21) and neglecting the second order terms, we obtain the 
following characterization of h'^ 



h'^ 



"=^t.H{2R\sm<Pp\)y'^^h'i-h'^] 

1 ^ 

^ N S [Gi( V2)(/^, - e'^^hi) + G2{cl>p/2)(hi - e'^^Vj)] . (22) 



=1 



In order to study the behavior of the perturbations /ij, we reduce the complexity of the 
problem, assuming that hj satisfles the following relation 



-imOn 



hj = C+(t)e*"^'^ + C-(t)e-*'"'^ h'^ = City""'' + C{t)e- 
Therefore, we can express hi in terms of hj as 



m G N. 
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Inserting the previous expressions in (22) and gathering terms in e -^^ and e 
characterize ^+ and ^_ as 



^^i^, we can 






where /i and /2 are defined in (14) and ( |15[ ). Through a simple manipulation of the sum for 
the linearized Cucker-Smale term, we obtain that the expression 



l^J]H{2R\smcl>,\)[e^^^ 



(m+l) 



^ 2 ^(2^1 ^i^^^l) [cos(0^(m ± 1)) - 1] + i^ J] iy(2i?| sin0^|) sin((/)^(m ± 1)), 

is real. Actually, H{2R\s\ii(j)p\) and sin((/)p(77i ± 1)) are respectively symmetric and anti- 
symmetric with respect to the values of 0^, so the imaginary part vanishes. Recalling the 
definition of 0^, we conclude 



1 ^ / 



. f2TTp 
in 



sm 



V N 



N/2 



N 



J]H(2R 



k^i 



sm 



/27rp 



( 2np , 

cosi — (m + l) ) -1 



)[sin2(^(m±l)) 



Therefore, we reduce the stability analysis to the following system 
^^'i\ /him) him) \ /e+\ /J+im) 



.r: 



yhirn) Ii{-m)l \^_ 







J-{m)l \t 



M 



Taking the conjugate in the second equation and relabeling ^_ with ^_ as in ^, the previous 
system is equivalent to 

1 

h{m) him) J+i^ 
\him) hi-m) J-im)j 



d_ 






Id 

M J 






(23) 



V+ 

\v-/ 

where 77+ = ^V- 

At this point, we will do a stability analysis based on the eigenvalues of the matrix of the 
previous system in a similar way as in Section 3.1 If instead of the self-propelled/ friction 
term we use a Cucker-Smale type alignment term 

N 



J]9{\xj-xi\){vj-vi), 



1^1 
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where g{r) denotes any strictly positive function, the corresponding stabihty matrix Lqs for 
the flock reads 

Id 



^^s - VM -G 

As before, M denotes stability matrix of the first order model. As the alignment term is linear 
in the velocity, the matrix G acts on v = {vi, . . . , vn)^ , '^j ^ I^^, according to the relation 

N 



In particular, if we denote \\vj — vi\\2 •= {^j — ^lYi^j ~ ^l) then this relation implies that 

N 

2 



1 ^ 
'G'v = - J] 9i\xj-xi\)\\vj-vi\\l. 



Consequently, G is positive semi-definite and Gv = if and only if v is "constant" in the sense 
that Vj = w for some fixed w e M?. In other words, ker(G) = span {61,62}. By translation 
invariance of the first order model, both 61 e ker(M) and 62 e ker(M) as well. 

Note that the eigenvalue problem for Lqs is again equivalent to the following quadratic 
eigenvalue problem for x e C^^: A^x + AGx — Mx = 0. Assuming the normalization x*x = 1, 
the quadratic formula then implies that 

_ -x*Gx ± V(x*Gx)2 + 4x*Mx 
^- 2 • 

From this relation and the fact that ker(G) c ker(M) we conclude that ker(Lcs) = ker(M), 
and moreover that 3ft(A) =0 <^ A = 0. Furthermore, ei and 62 generate a single gener- 
alized eigenvector whereas each remaining x e ker(M) generates no generalized eigenvectors. 
Indeed, corresponding to each x e ker(M) the system of equations 

Id \ /u\ ^ /x 
M -GJ [wj ~ \0 

has a solution if and only if Gx = and u e ker(M). Additionally, if x = e^ then for any 
u G ker(M) the system of equations 

Id \ /u\ ^ /u 
M -GJ [wj ~ [y: 

has no solutions. This follows by multiplying the second equation by x*, then using the facts 
that 6^ G ker(G) ^ ker(M) and the facts that G and M are symmetric. In other words, if 
g{r) is any strictly positive function then 

det(Lcs - Aid) = A2+^i^(ker(M))^^^^)^ 

for some polynomial Pg{X) that has non-zero roots. Since this equation holds for any strictly 
positive function ^(r), we may follow the proof of Theorem [l] to conclude that the second 
order model has an eigenvalue with positive real part if and only if the first order system 
has a positive eigenvalue. Moreover, the vectors (6^,6^) for each z = 1,2 furnish generalized 
eigenvectors with eigenvalue zero, so the ring flock is always linearly unstable. As a summary, 
we have shown: 
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Theorem 2. The linearized second order system ^ around the flock ring solution has an 
eigenvalue with positive real part if and only if the linearized first order system around the 
ring solution has a positive eigenvalue. As a consequence, the flock ring solution is unstable 
for m-mode perturbations for the second order model ^ if and only if the ring solution is 
unstable for m-mode perturbations for the first order model Q. 

The hnear stabihty analysis of the previous system leads to the characterization of the 
stability areas in Figure [6) We show that the stability parameter regions for different values 
of N^ and 7 = 1 coincide with the ones in Figure [3) 




Figure 6. Stability regions for flock ring solutions with the Cucker-Smale 
term for different values of N. 



We also investigate the behavior of the eigenvalue with the largest real part, 3?(Ai), of 
the linearized system (23) against the increasing value of communication strength 7. In 



Figure [7| as the potential gets more repulsive at the origin, we see the change from stability 
to instability, and the rate of convergence to equilibrium depending on 7. 




Figure 7. The magnitude of 5ft(Ai) is influenced by 7, for different values of 
b and fixed a = 5, N = 10000. 



4. Stability for mill solutions 

This section is devoted to complement the results in [6J by analyzing the stability of mill 
ring solutions with repulsion. 
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4.1. Linear stability analysis. Let us consider the transformation 



U{t) = Oit)xjit) 
\zjit) = Oit)vjit) 

where 0{t) is the rotation matrix defined as 

-U 



,A^ 



0{t) 



s = 



and d{t) = Se 



St 



Evaluating yj{t) and Zj{t) and after some straightforward computations, we get 

^ ijjit) = Syj{t) + Zj{t) 
\ zj{t) = Szj{t) + ia-(3\zj\')zjit) - J]vWiyi-v.) > J = h...,N. 



N 



Vj) 



1=1 

1^3 



A hnear stabihty analysis for mill rings was performed in [6J . Actually, for a fixed number of 
particles, we have a mill ring solution given by {y^j,z^A = (i?e*^-', 0), where Oj = -^, and R 
determined by equation ^ . With the same notation as in Section |3J the analysis in [6J leads 
to the linear system 

AV\ / 








-ujia + c<j + 



him) 



\ 



bjia 



him) 



ujia 



-ujia 



him) 
Mi(-m) 



1 


-a — 2uji 
a 



\ fcA 


1 


C- 


a 


V+ 


a + 2ujiJ \ri-J 



(24) 



Let us remind that the perturbations are of the form yj{t) = Re^^^{l + hj{t)), with hj 

2, 3, ... , such that \hj\ « 1 and satisfying ([T2|), with (77+, 77_^ 



-imOn 



m 



(C+^C-)- We will make use of (24) to study the stability of mill rings with repulsion. 



4.2. Numerical tests. Unlike the case of flock solutions where the asymptotic speed does 
not play role in the linear stability we will show that the asymptotic speed \uo\ can be used 
as a bifurcation parameter for mills. 

In Table [3] we numerically investigate the behavior of the stability region for a fixed number 
of particles, N = 1000, and for increasing values of the asymptotic speed \uq\. We observe 
that the stability region shrinks with respect to a and gets larger with respect to b. Each 
stability region in Table [3] is computed out of the intersection of the stable areas for the 
system (24) for each perturbation mode m ^ 2. Note that for |^xo| = the stability region 
coincides with the one for the the first order model ^ and for the flock ring solution. 

A similar analysis, as done in Subsection 3^, can be performed to study the formation of 
fat mills and clustered mill solutions. We show how both the fattening and the clustering 
instability are triggered by tunning the asymptotic speed for a choice of the interaction 
potential (a and b). 

In the case of flock ring solutions we observe cluster solutions or annulus solutions when 
parameters a and b are chosen respectively "below" or "above" the stability region. In the 
case of mill solutions, a similar behavior is observed, but this will depend also on the chosen 
value of l^xol- As an example, we fix (a, b) = (5, 0.5), marked as (0) in Table|3J and we observe 
the behavioral change of the system for increasing values of the asymptotic speed. 

Table [4] exhibits this switching behavior from a fat mill to a cluster pattern along with the 
increment of the asymptotic speed. We observe that for small values of the asymptotic speed 
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Table 3. Stability region for N = 1000 and different values of the asymptotic 
speed I^xqI- Markers (0) and (^) correspond to the explored parameters in 
Table 11 and Tabled 



fat mill solutions are stable patterns, but when increasing the value of |i/o| the stable solutions 
form a clustered mill. The first row shows the evolution of the system with asymptotic speed 
\uq\ = 0.25 towards an annulus mill. In the second row we take |^xo| = 0.5 the previous 
stable solution is reshaped to a fat clusters pattern. The speed in the third row is switched 
to I^xqI = 5 and clusters on lines emerge as a stable configuration. Increasing the speed to 
l^xol = 50 in the fourth row we can observe that clusters on "points" are stable solutions. 

In Figure [s] we numerically show how the stability region looks like in terms of {\uo\^b)^ 
with a fixed at 0.5, and we enlighten with marker points (0) the parameter choices of Table[4| 
first and second lines. 

For the sake of completeness, we enrich the analysis fixing \uq\ = 0.5 and considering 
different values of 6, in order to cross the stability region. Therefore in Table [5] we show the 
evolution of a mill ring solution with b taken subsequently equal to 0.5, 1.25,3.5, parameter 
choices are marked as {^) in Figure |8J The first line of Table [5] shows the convergence to the 
same stable state as the one in second line of Table [4| but, since the system starts to evolve 
directly from a ring mill solution, the transient behavior is different. Parameters in second 
line belong to the stability region, see Figure [8) Therefore, the stable state becomes a mill 
ring solution. Finally, in the third line we increase b and a three point cluster solution is 
observed as stable pattern. 



4.2.1. Mill to Flock and Flock to Mill behavior. We numerically investigate the stability of 
mill and flock ring solutions for small values of the asymptotic speed, \uo\, and the parameter 
6, which corresponds to a strong repulsion condition. 



\uo\ 



0.25 



0.5 



50 
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Stable 





m^ 



f 






^ 



) \ 

Table 4. A^ = 1000 particles, a = 5,6 = 0.5. The table shows the evolution 
of a mill ring for increasing values of the speed |i/o|- Each row depicts the 
behavior of the system for a fixed speed, until a stable state is reached. The 
evolution of the second, third and fourth row is computed starting from the 
stable pattern of the previous line. 
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Figure 8. Stability region for mill ring solution with parameter a = 5, 
N = 1000. Markers (0) and {^) depict the parameter choices respectively 
for Table [land Table [5l 
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Table 5. A^ = 1000 particles, a = 5, |ixo| = 0.5. The table shows the evolution 
of a mill ring for increasing values of 6, i.e. decreasing repulsion. The evolution 
of the second and the third row is computed starting from the stable pattern 
of the previous line. 



We perform two representative simulations showing that for a particular choice of the 
parameters, mill ring solutions can switch to fat flock solutions and conversely flock mill 
solutions switch to fat mill patterns. 
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Table 6. System with A^ = 100 agents, parameters are fixed a = 4 and 
h = 0.0005 and \u^\ = 0.01. The first row shows that the initial mill ring 
configuration is unstable. The second row outlines the self organization of the 
system in a fat flock configuration. 



In Table |6] we take N = 100 particles and we fix a = 4,6 = 0.0005 and \uo\ = 0.01. The 
frames in the first row show the instability of mill ring solutions for this choice of parameters. 
The system initially evolves to an almost chaotic state, then particles start to organize rotating 
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around the center of mass. This rotation actually causes the alignment of the agents and the 
final fat flock configuration described in the second row. In Table [7| we consider as initial 
state a flock ring solution. The parameters of the model are N = 100, a = 4, b = 0.001 
and \uo\ = 0.1. The first row of the table illustrates that the initial configuration is not a 
stable solution. Therefore, the symmetry of the flock ring is broken and the system exhibits 
a chaotic behavior. In the second row a rotating dynamic emerges out of the disordered state 
and finally the system stabilizes to a fat mill solution. 






• i~ 
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Table 7. System with A^ = 100 agents, parameters a = 4 and h = 0.001 and 
l^xol = 0.1. The first row shows the instability of the flock ring solution while 
the second exhibits the convergence to a fat mill type solution. 



These numerical tests show surprisingly that it is possible, with a particular choice of the 
parameters, to obtain mill configurations out of perturbations of initial flock solutions and 
flock solutions out of perturbations of mill ring solutions. These heteroclinic-kind solutions 
have not been previously reported. We also remark that the parameter choice is connected to 
the number of agents we are considering; changing N means finding another set of parameters 
for which the same switching behavior occurs. 
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